Poster project

Data: Load data

library(st514)
# Indlæsning af data
data("T6-19")
T6.19 <- tbl
T6.19
##       V1    V2 V3
## 1  271.0 18.50  F
## 2  477.0 82.50  F
## 3  306.3 23.40  F
## 4  365.3 33.50  F
## 5  466.0 69.00  F
## 6  440.7 54.00  F
## 7  315.0 24.97  F
## 8  417.5 56.75  F
## 9  307.3 23.15  F
## 10 319.0 29.51  F
## 11 303.9 19.98  F
## 12 331.7 24.00  F
## 13 435.0 70.37  F
## 14 261.3 15.50  F
## 15 384.8 63.00  F
## 16 360.3 39.00  F
## 17 441.4 53.00  F
## 18 246.7 15.75  F
## 19 365.3 44.00  F
## 20 336.8 30.00  F
## 21 326.7 34.00  F
## 22 312.0 25.00  F
## 23 226.7  9.25  F
## 24 347.4 30.00  F
## 25 280.2 15.25  F
## 26 290.7 21.50  F
## 27 438.6 57.00  F
## 28 377.1 61.50  F
## 29 176.7  3.00  M
## 30 259.5  9.75  M
## 31 258.0 10.07  M
## 32 229.8  7.50  M
## 33 233.0  6.25  M
## 34 237.5  9.85  M
## 35 268.3 10.00  M
## 36 222.5  9.00  M
## 37 186.5  3.75  M
## 38 238.8  9.75  M
## 39 257.6  9.75  M
## 40 172.0  3.00  M
## 41 244.7 10.00  M
## 42 224.7  7.25  M
## 43 231.7  9.25  M
## 44 235.9  7.50  M
## 45 236.5  5.75  M
## 46 247.4  7.75  M
## 47 223.0  5.75  M
## 48 223.7  5.75  M
## 49 212.5  7.65  M
## 50 223.2  7.75  M
## 51 225.0  5.84  M
## 52 228.0  7.53  M
## 53 215.6  5.75  M
## 54 221.0  6.45  M
## 55 236.7  6.49  M
## 56 235.3  6.00  M

Data: update factorial strings

colnames(T6.19) <- c("Snout Vent Length (cm)", "Weight (kg)", "Gender")
levels(T6.19$Gender) <- c("Female", "Male")
T6.19
##    Snout Vent Length (cm) Weight (kg) Gender
## 1                   271.0       18.50 Female
## 2                   477.0       82.50 Female
## 3                   306.3       23.40 Female
## 4                   365.3       33.50 Female
## 5                   466.0       69.00 Female
## 6                   440.7       54.00 Female
## 7                   315.0       24.97 Female
## 8                   417.5       56.75 Female
## 9                   307.3       23.15 Female
## 10                  319.0       29.51 Female
## 11                  303.9       19.98 Female
## 12                  331.7       24.00 Female
## 13                  435.0       70.37 Female
## 14                  261.3       15.50 Female
## 15                  384.8       63.00 Female
## 16                  360.3       39.00 Female
## 17                  441.4       53.00 Female
## 18                  246.7       15.75 Female
## 19                  365.3       44.00 Female
## 20                  336.8       30.00 Female
## 21                  326.7       34.00 Female
## 22                  312.0       25.00 Female
## 23                  226.7        9.25 Female
## 24                  347.4       30.00 Female
## 25                  280.2       15.25 Female
## 26                  290.7       21.50 Female
## 27                  438.6       57.00 Female
## 28                  377.1       61.50 Female
## 29                  176.7        3.00   Male
## 30                  259.5        9.75   Male
## 31                  258.0       10.07   Male
## 32                  229.8        7.50   Male
## 33                  233.0        6.25   Male
## 34                  237.5        9.85   Male
## 35                  268.3       10.00   Male
## 36                  222.5        9.00   Male
## 37                  186.5        3.75   Male
## 38                  238.8        9.75   Male
## 39                  257.6        9.75   Male
## 40                  172.0        3.00   Male
## 41                  244.7       10.00   Male
## 42                  224.7        7.25   Male
## 43                  231.7        9.25   Male
## 44                  235.9        7.50   Male
## 45                  236.5        5.75   Male
## 46                  247.4        7.75   Male
## 47                  223.0        5.75   Male
## 48                  223.7        5.75   Male
## 49                  212.5        7.65   Male
## 50                  223.2        7.75   Male
## 51                  225.0        5.84   Male
## 52                  228.0        7.53   Male
## 53                  215.6        5.75   Male
## 54                  221.0        6.45   Male
## 55                  236.7        6.49   Male
## 56                  235.3        6.00   Male

Data: Working with the data

mean(T6.19$"Snout Vent Length (cm)")
## [1] 288.5143
mean(T6.19$"Weight (kg)")
## [1] 22.27696
# Variance of x1 (using n divisor)
var(T6.19$"Snout Vent Length (cm)") * (nrow(T6.19) - 1) / nrow(T6.19)
## [1] 6084.646
# Variance of x2 (using n divisor)
var(T6.19$"Weight (kg)") * (nrow(T6.19) - 1) / nrow(T6.19)
## [1] 421.5405
# Covariance between x1 and x2 (using n divisor)
cov(T6.19$"Snout Vent Length (cm)", T6.19$"Weight (kg)") * (nrow(T6.19) - 1) / nrow(T6.19)
## [1] 1539.699
# Correlation between x1 and x2
cor(T6.19$"Snout Vent Length (cm)", T6.19$"Weight (kg)")
## [1] 0.9613875
# Create a scatterplot with a marginal dot plot for the variables x1 and x2
marginal.dot.plot(T6.19$"Snout Vent Length (cm)", T6.19$"Weight (kg)", xlab = expression("Snout Vent Length (cm) x"[1]), ylab = expression("Weight (kg) x"[2]))

# Variance-covariance matrix:

#create t6.19 nuerical withouht column gender
T6.19_numerical <- T6.19[, 1:2]

#print(round(cov(T6.19_numerical) * (nrow(T6.19) - 1) / nrow(T6.19)),3)
print(round(cov(T6.19_numerical) * (nrow(T6.19) - 1) / nrow(T6.19), digits=3))
##                        Snout Vent Length (cm) Weight (kg)
## Snout Vent Length (cm)               6084.646    1539.699
## Weight (kg)                          1539.699     421.540
# Find standardafvigelser (kvadratrod af diagonale værdier i varians-kovarians-matrixen)
vcov_matrix <- cov(T6.19_numerical) * (nrow(T6.19) - 1) / nrow(T6.19)

# Træk diagonalen ud (varians for hver variabel), og tag kvadratroden
std_devs <- sqrt(diag(vcov_matrix))

# Udskriv med beskrivende tekst
cat("Standardafvigelser (√varians):\n")
## Standardafvigelser (√varians):
print(round(std_devs, 2))
## Snout Vent Length (cm)            Weight (kg) 
##                  78.00                  20.53

DATA: Explanatory Data Analysis

# Dimensioner
dim(T6.19)
## [1] 56  3

56 observations and 3 variables

#Variabeltyper
str(T6.19)
## 'data.frame':    56 obs. of  3 variables:
##  $ Snout Vent Length (cm): num  271 477 306 365 466 ...
##  $ Weight (kg)           : num  18.5 82.5 23.4 33.5 69 ...
##  $ Gender                : Factor w/ 2 levels "Female","Male": 1 1 1 1 1 1 1 1 1 1 ...
#Sammendrag
summary(T6.19)
##  Snout Vent Length (cm)  Weight (kg)       Gender  
##  Min.   :172.0          Min.   : 3.00   Female:28  
##  1st Qu.:229.3          1st Qu.: 7.50   Male  :28  
##  Median :258.8          Median :10.04              
##  Mean   :288.5          Mean   :22.28              
##  3rd Qu.:333.0          3rd Qu.:30.00              
##  Max.   :477.0          Max.   :82.50

EDA: Søjlediagram

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(ggplot2)

# Tilføj et unikt ID til hver observation
T6.19_obs <- T6.19 %>%
  mutate(Obs = paste0("Obs", row_number()))

# Gør dataen lang for Snout og Weight
T6.19_long <- T6.19_obs %>%
  dplyr::select(`Obs`, `Gender`, `Snout Vent Length (cm)`, `Weight (kg)`) %>%
  pivot_longer(cols = c(`Snout Vent Length (cm)`, `Weight (kg)`),
               names_to = "Variable", values_to = "Value")

# Standardiser variabelnavne lidt for pænere labels
library(stringr)

T6.19_long$Variable <- str_replace_all(T6.19_long$Variable,
                                       c("Snout Vent Length \\(cm\\)" = "Snout Length",
                                         "Weight \\(kg\\)" = "Weight"))

# Plot: Hver søjle er én observation
ggplot(T6.19_long, aes(x = Obs, y = Value, fill = Gender)) +
  geom_bar(stat = "identity") +
  facet_wrap(~ Variable, scales = "free_y") +
  labs(
    title = "Snout Length og Weight for alle observationer",
    x = "Observationer", y = "Værdi"
  ) +
  scale_fill_manual(values = c("Female" = "#F8766D", "Male" = "#619CFF")) +
  theme_minimal() +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

library(dplyr)
library(tidyr)
library(ggplot2)
library(stringr)

# Tilføj numerisk ID og behold som karakter for rækkefølge
T6.19_obs <- T6.19 %>%
  mutate(Obs = row_number())

# Gør data lang
T6.19_long <- T6.19_obs %>%
  dplyr::select(Obs, Gender, `Snout Vent Length (cm)`, `Weight (kg)`) %>%
  pivot_longer(cols = c(`Snout Vent Length (cm)`, `Weight (kg)`),
               names_to = "Variable", values_to = "Value")

# Standardiser navne
T6.19_long$Variable <- str_replace_all(T6.19_long$Variable,
                                       c("Snout Vent Length \\(cm\\)" = "Snout Length",
                                         "Weight \\(kg\\)" = "Weight"))

# Gør Obs til faktor for rækkefølge
T6.19_long$Obs <- factor(T6.19_long$Obs, levels = 1:56)

# Kombiner observation og variabel til unik bar (bruges ikke i plot her men gemmes til evt. brug)
T6.19_long$ObsVar <- interaction(T6.19_long$Obs, T6.19_long$Variable, sep = "_")

# Plot begge variabler som separate søjler i ét diagram
ggplot(T6.19_long, aes(x = Obs, y = Value, fill = Gender)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9), aes(group = Variable)) +
  labs(
    title = "Snout Length og Weight pr. observation",
    x = "Observation", y = "Værdi"
  ) +
  scale_fill_manual(values = c("Female" = "#F8766D", "Male" = "#619CFF")) +
  theme_minimal()

library(dplyr)
library(tidyr)
library(ggplot2)

# Beregn gennemsnit for hver kombination af variabel og køn
T6.19_means <- T6.19 %>%
  group_by(Gender) %>%
  summarise(
    Snout = mean(`Snout Vent Length (cm)`),
    Weight = mean(`Weight (kg)`)
  ) %>%
  pivot_longer(cols = c(Snout, Weight), names_to = "Variable", values_to = "Mean")

# Plot
ggplot(T6.19_means, aes(x = Gender, y = Mean, fill = Variable)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(
    title = "Gennemsnitlig Snout Length og Weight fordelt på køn",
    x = "Køn", y = "Gennemsnitlig værdi"
  ) +
  scale_fill_manual(values = c("Snout" = "#619CFF", "Weight" = "#F8766D")) +
  theme_minimal()

# Scatterplot of length and weight given gender=female

## Scatterplot: Females only
T6.19_female <- subset(T6.19, Gender == "Female")

plot(T6.19_female$`Snout Vent Length (cm)`, T6.19_female$`Weight (kg)`,
     xlab = expression("Snout Vent Length (cm) x"[1]),
     ylab = expression("Weight (kg) x"[2]),
     main = "Scatterplot of Female Anacondas",
     pch = 19, col = "darkred")

## Scatterplot: Males only
T6.19_male <- subset(T6.19, Gender == "Male")

plot(T6.19_male$`Snout Vent Length (cm)`, T6.19_male$`Weight (kg)`,
     xlab = expression("Snout Vent Length (cm) x"[1]),
     ylab = expression("Weight (kg) x"[2]),
     main = "Scatterplot of Male Anacondas",
     pch = 19, col = "darkblue")

plot(T6.19$`Snout Vent Length (cm)`, T6.19$`Weight (kg)`,
     xlab = expression("Snout Vent Length (cm) x"[1]),
     ylab = expression("Weight (kg) x"[2]),
     main = "Scatterplot by Gender",
     pch = 19,
     col = ifelse(T6.19$Gender == "Female", "darkred", "darkblue"))

legend("topleft", legend = c("Female", "Male"),
       col = c("darkred", "darkblue"), pch = 19)

#Mean vector
print(colMeans(T6.19_numerical, na.rm = TRUE), digits = 5)
## Snout Vent Length (cm)            Weight (kg) 
##                288.514                 22.277
# colMeans of T6.19 where gender=female
# Column means for females
T6.19_female <- subset(T6.19, Gender == "Female")
T6.19_female_numerical <- T6.19_female[, 1:2]
print(colMeans(T6.19_female_numerical, na.rm = TRUE), digits = 5)
## Snout Vent Length (cm)            Weight (kg) 
##                348.275                 37.264
# colMeans of T6.19 where gender=male
# Column means for males
T6.19_male <- subset(T6.19, Gender == "Male")
T6.19_male_numerical <- T6.19_male[, 1:2]
print(colMeans(T6.19_male_numerical, na.rm = TRUE), digits = 5)
## Snout Vent Length (cm)            Weight (kg) 
##               228.7536                 7.2904
#correlation matrix
print(cor(T6.19_numerical))
##                        Snout Vent Length (cm) Weight (kg)
## Snout Vent Length (cm)              1.0000000   0.9613875
## Weight (kg)                         0.9613875   1.0000000

Pairsplot

pairs(T6.19[, 1:2], 
      main = "Scatterplot Matrix: Length and Weight",
      pch = 21,
      bg = c("lightpink", "lightblue")[T6.19$Gender])

#pairs(T6.19_numerical, diag.panel = panel.boxplot, labels = "")
pairs(T6.19_numerical,
      diag.panel = panel.boxplot,
      labels = c("Snout Length", "Weight"),
      main = "Scatterplot Matrix: Length and Weight")
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter

# Calculate correlation coefficient
cat("Correlation coefficient:", cor(T6.19$"Snout Vent Length (cm)", T6.19$"Weight (kg)"), "\n")
## Correlation coefficient: 0.9613875
# Calculate correlation coefficient
cat("Correlation coefficient:", cor(T6.19$"Weight (kg)", T6.19$"Snout Vent Length (cm)"), "\n")
## Correlation coefficient: 0.9613875
# Plot af data
plot(T6.19_numerical,
     xlim = c(0, 500),
     ylim = c(0, 100),
     main = "Scatterplot med stikprøvegennemsnit",
     xlab = "Snout Vent Length (cm) x1", ylab = "Weight (kg) x2", pch = 16)

grid()


# Beregning af stikprøvegennemsnit
x_bar <- colMeans(T6.19_numerical)
#x_bar


# Tilføj stikprøvegennemsnit til plot
points(x_bar[1], x_bar[2], pch = "X", col = "red", cex = 2)

# Tilføj en forklaring
legend("topleft", 
       legend = c("56 Observationer", "Stikprøvegennemsnit"),
       pch = c(16, 4), 
       col = c("black", "red"),
       pt.lwd = c(1, 3))

# Vis stikprøvegennemsnittet
cat("Stikprøvegennemsnittet er:", x_bar, "\n")
## Stikprøvegennemsnittet er: 288.5143 22.27696
# Udtræk data
Xf <- T6.19[T6.19$Gender == "Female", 1:2]
Xm <- T6.19[T6.19$Gender == "Male", 1:2]
n1 <- nrow(Xf); n2 <- nrow(Xm)
xbar1 <- colMeans(Xf)
xbar2 <- colMeans(Xm)
diff_mean <- xbar1 - xbar2

# Pooled kovariansmatrix
S1 <- cov(Xf); S2 <- cov(Xm)
Spooled <- ((n1 - 1)*S1 + (n2 - 1)*S2)/(n1 + n2 - 2)

# Egenværdier og egenvektorer
eig <- eigen(Spooled)
lambda <- eig$values      # egenværdier
vectors <- eig$vectors    # egenvektorer

# Parametre for ellipse
alpha <- 0.05
p <- 2
n <- (n1 * n2) / (n1 + n2)
T2_crit <- ((n1 + n2 - 2) * p / (n1 + n2 - p - 1)) * qf(1 - alpha, p, n1 + n2 - p - 1)

# Generér ellipsepunkter
theta <- seq(0, 2*pi, length = 100)
circle <- rbind(cos(theta), sin(theta))
scaling <- sqrt(T2_crit / n) * sqrt(lambda)
ellipse_points <- t(vectors %*% diag(scaling) %*% circle + diff_mean)

# Plot
plot(ellipse_points, type = "l", asp = 1, lwd = 2,
     xlab = "Snout Vent Length (cm)", ylab = "Weight (kg)",
     main = "Hotelling’s T² Ellipse (uden ellipse-pakke)")
points(0, 0, pch = 4, col = "red")
text(0, 0, "0", pos = 1, cex = 0.8)
points(diff_mean[1], diff_mean[2], pch = 19, col = "blue")
text(diff_mean[1], diff_mean[2], "M1 - M2", pos = 4, cex = 0.8)

# Beregn generaliseret stikprøvevarians for Spooled matrix
Xf <- T6.19[T6.19$Gender == "Female", 1:2]
Xm <- T6.19[T6.19$Gender == "Male", 1:2]
n1 <- nrow(Xf); n2 <- nrow(Xm)
S1 <- cov(Xf); S2 <- cov(Xm)
Spooled <- ((n1 - 1)*S1 + (n2 - 1)*S2)/(n1 + n2 - 2)

# Determinant = generaliseret varians
gen_var <- det(Spooled)
cat("Generaliseret stikprøvevarians (|S|):", round(gen_var, 2))
## Generaliseret stikprøvevarians (|S|): 86170.38

Test for normalitet (før hypotesetests)

1.1 Univariate normalitetstest (marginal):

Brug Q–Q plot for hver af de to variable (Snout Vent Length og Weight).

F.eks. som i Øvelse 5.4.B.

1.2 Multivariat normalitetstest:

Brug Q–Q plot for Mahalanobis afstande mod kvantilfunktion for χ² (chi-square plot).

F.eks. baseret på procedure i bogens Kapitel 4, Ex. 4.13 og Resultat 4.7
# Lav Q-Q plots for hver af de tre variable
par(mfrow = c(2, 2))

for (i in 1:p) {
  qqnorm(T6.19_numerical[,i], main = colnames(T6.19_numerical)[i])
  qqline(T6.19_numerical[,i])
}

# Beregn Mahalanobis-afstandene
Sx <- cov(T6.19_numerical)
D2 <- mahalanobis(T6.19_numerical, colMeans(T6.19_numerical), Sx)

# Lav et chi-i-anden Q-Q plot for multivariate normalitet
qqplot(qchisq(ppoints(n, a = 0.5), df = p), D2,
       ylab = "Mahalanobis afstande",
       xlab = bquote("Kvantiler af " ~ chi[.(p)]^2),
       main = bquote("Q-Q plot af Mahalanobis" * ~ D^2 * 
                       " vs. kvantiler af" * ~ chi[.(p)]^2))
abline(0, 1)

library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
library(car) # til powerTransform
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
# Antal observationer og dimensioner
n <- nrow(T6.19_numerical)
p <- ncol(T6.19_numerical)

# Box-Cox transformationer
pt_result <- powerTransform(T6.19_numerical)
lambda_vec <- pt_result$lambda
print(lambda_vec)
## Snout Vent Length (cm)            Weight (kg) 
##            -0.49978201            -0.03326786
# Transformer hver variabel
T6.19_transformed <- as.data.frame(
  sapply(1:p, function(i) {
    lambda <- lambda_vec[i]
    x <- T6.19_numerical[, i]
    if (lambda == 0) log(x) else (x^lambda - 1) / lambda
  })
)
colnames(T6.19_transformed) <- colnames(T6.19_numerical)

# Plot QQ-plots for de transformerede variable
par(mfrow = c(2, 2))
for (i in 1:p) {
  qqnorm(T6.19_transformed[, i], 
         main = paste("Q-Q plot (transformeret):", colnames(T6.19_transformed)[i]))
  qqline(T6.19_transformed[, i])
}

# Mahalanobis-afstande (baseret på transformerede data)
Sx <- cov(T6.19_transformed)
D2 <- mahalanobis(T6.19_transformed, colMeans(T6.19_transformed), Sx)

# Q-Q plot af Mahalanobis afstande vs. chi^2
qqplot(qchisq(ppoints(n, a = 0.5), df = p), D2,
       ylab = expression("Mahalanobis afstande"),
       xlab = bquote("Kvantiler af " ~ chi[.(p)]^2),
       main = bquote("Q-Q plot af Mahalanobis" * ~ D^2 * " (transformeret)"))
abline(0, 1)

.

2. Univariate inferens (T-tests)

2. Univariate analyser

  • T-test for:
# Univariate t-tests for Snout Vent Length og Weight mod Gender

# Antag at køn er korrekt kodet som faktor med levels: Female, Male
T6.19$Gender <- factor(T6.19$Gender, levels = c("Female", "Male"))
  • længde vs. køn
# T-test for Snout Vent Length (cm) mellem køn
t_snout <- t.test(`Snout Vent Length (cm)` ~ Gender, data = T6.19, var.equal = TRUE)

#var.equal = TRUE bruges her i tråd med forudsætningen fra tidligere øvelser, hvor vi antager lig varians ved sammenligninger af grupper.


print(t_snout)
## 
##  Two Sample t-test
## 
## data:  Snout Vent Length (cm) by Gender
## t = 8.7597, df = 54, p-value = 0.000000000005976
## alternative hypothesis: true difference in means between group Female and group Male is not equal to 0
## 95 percent confidence interval:
##   92.16599 146.87687
## sample estimates:
## mean in group Female   mean in group Male 
##             348.2750             228.7536

✅ Korrekt fortolkning:

H₀ (nulhypotese): Der er ingen forskel i gennemsnitlig snout length mellem han- og hun-anacondaer.

H₁ (alternativ hypotese): Der er en forskel i gennemsnitlig snout length mellem kønnene.

📌 Resultat:

P-værdien er ekstremt lille (≈ 5.98e-12), dvs. langt under den typiske signifikansgrænse (α = 0.05).

Derfor forkaster vi H₀: Der er statistisk evidens for, at snout length afhænger af køn.

Konfidensintervallet [92.17, 146.88][92.17, 146.88] indeholder ikke nul, hvilket styrker denne konklusion.

Hun-anacondaer (♀) har i gennemsnit ca. 120 cm længere snout end hanner (♂).

❗ Korrigeret fortolkning:

"Da p-værdien er meget lav (<< 0.05), forkastes H₀. Der er altså en statistisk signifikant forskel i snout length mellem hun- og han-anacondaer. Det tyder på, at snout length afhænger af køn."

💥 OBS: Din formulering “H0 hypotesen ikke forkastes baseret på lav p-værdi” er forkert — lav p-værdi betyder netop, at vi forkaster H₀.

  • vægt vs. køn
# T-test for Weight (kg) mellem køn
t_weight <- t.test(`Weight (kg)` ~ Gender, data = T6.19, var.equal = TRUE)

#var.equal = TRUE bruges her i tråd med forudsætningen fra tidligere øvelser, hvor vi antager lig varians ved sammenligninger af grupper.


print(t_weight)
## 
##  Two Sample t-test
## 
## data:  Weight (kg) by Gender
## t = 7.8475, df = 54, p-value = 0.0000000001739
## alternative hypothesis: true difference in means between group Female and group Male is not equal to 0
## 95 percent confidence interval:
##  22.31565 37.63078
## sample estimates:
## mean in group Female   mean in group Male 
##            37.263571             7.290357

✅ Korrekt fortolkning:

H₀ (nulhypotese): Der er ingen forskel i gennemsnitlig vægt mellem kønnene.

H₁ (alternativ hypotese): Der er forskel i gennemsnitlig vægt mellem kønnene.

📌 Resultat:

P-værdien er ekstremt lav (≈ 1.74e-10), altså klart under 0.05 → forkast H₀.

Der er statistisk belæg for, at køn har stor betydning for vægt.

Konfidensintervallet viser, at hunner vejer 22–38 kg mere end hanner i gennemsnit.

❗ Korrigeret fortolkning:

"P-værdien er meget lav, hvilket betyder, at vi forkaster H₀. Der er statistisk signifikant forskel i vægt mellem køn. Vægten afhænger af køn, og hun-anacondaer vejer gennemsnitligt betydeligt mere."

(Sammenfatning)

✅ Samlet konklusion:

Test p-værdi H₀ forkastes? Fortolkning Snout Length ~ Gender 5.98e-12 Ja Snout length afhænger af køn Weight ~ Gender 1.74e-10 Ja Vægt afhænger af køn

📊 Tabel til posteren (Univariate T-tests)

Variabel Middelværdi (♀) Middelværdi (♂) T-værdi P-værdi 95% CI for forskel Konklusion
Snout Length (cm) 348.28 228.75 8.76 < 0.000000000006 [92.17, 146.88] Signifikant forskel. Afhænger af køn.
Weight (kg) 37.26 7.29 7.85 < 0.000000000174 [22.32, 37.63] Signifikant forskel. Afhænger af køn.

🧾 Noter til posteren (fortolkning og forklaring):

Univariate T-tests blev udført for hver af de to variable:

  Snout Length (cm) og Weight (kg) fordelt på køn (♀ vs ♂).

  Begge tests viser ekstremt lave p-værdier (<< 0.05), hvilket betyder, at vi forkaster nulhypotesen.

🟢 Fortolkning:

  Der er statistisk signifikant forskel i både længde og vægt mellem hun- og han-anacondaer.

  Hun-anacondaer har i gennemsnit 120 cm længere snout og 30 kg højere vægt end hanner.

  Disse forskelle er ikke blot statistisk signifikante, men også biologisk relevante.
  
  
  

📦 Boxplots for Snout Length og Weight fordelt på køn

library(ggplot2)
library(tidyr)
library(dplyr)

# Gør data lang
T6.19_long <- T6.19 %>%
  dplyr::select(Gender, `Snout Vent Length (cm)`, `Weight (kg)`) %>%
  pivot_longer(cols = -Gender, names_to = "Variable", values_to = "Value")

# Standardiser labels
T6.19_long$Variable <- gsub("Snout Vent Length \\(cm\\)", "Snout Length", T6.19_long$Variable)
T6.19_long$Variable <- gsub("Weight \\(kg\\)", "Weight", T6.19_long$Variable)

# Boxplot
ggplot(T6.19_long, aes(x = Gender, y = Value, fill = Gender)) +
  geom_boxplot(alpha = 0.7) +
  facet_wrap(~ Variable, scales = "free_y") +
  labs(
    title = "Boxplots af Snout Length og Weight fordelt på køn",
    x = "Køn", y = "Værdi"
  ) +
  scale_fill_manual(values = c("Female" = "#F8766D", "Male" = "#619CFF")) +
  theme_minimal()

# Violinplot
ggplot(T6.19_long, aes(x = Gender, y = Value, fill = Gender)) +
  geom_violin(trim = FALSE, alpha = 0.6, bw = 10) +  # Juster bw (bandwidth)
  geom_boxplot(width = 0.1, fill = "white", outlier.shape = NA) +
  facet_wrap(~ Variable, scales = "free_y") +
  labs(
    title = "Violinplots af Snout Length og Weight fordelt på køn",
    x = "Køn", y = "Værdi"
  ) +
  scale_fill_manual(values = c("Female" = "#F8766D", "Male" = "#619CFF")) +
  theme_minimal()

ggplot(T6.19_long, aes(x = Gender, y = Value, fill = Gender)) +
  geom_violin(trim = FALSE, scale = "count", alpha = 0.6) +
  geom_boxplot(width = 0.1, fill = "white", outlier.shape = NA) +
  facet_wrap(~ Variable, scales = "free_y") +
  labs(
    title = "Violinplots af Snout Length og Weight fordelt på køn",
    x = "Køn", y = "Værdi"
  ) +
  scale_fill_manual(values = c("Female" = "#F8766D", "Male" = "#619CFF")) +
  theme_minimal()

# Violinplot: Snout Length
ggplot(T6.19, aes(x = Gender, y = `Snout Vent Length (cm)`, fill = Gender)) +
  geom_violin(trim = FALSE, alpha = 0.6) +
  geom_boxplot(width = 0.1, fill = "white", outlier.shape = NA) +
  labs(
    title = "Violin Plot: Snout Length by Gender",
    x = "Gender", y = "Snout Length (cm)"
  ) +
  scale_fill_manual(values = c("Female" = "#F8766D", "Male" = "#619CFF")) +
  theme_minimal()

# Violinplot: Weight
ggplot(T6.19, aes(x = Gender, y = `Weight (kg)`, fill = Gender)) +
  geom_violin(trim = FALSE, alpha = 0.6) +
  geom_boxplot(width = 0.1, fill = "white", outlier.shape = NA) +
  labs(
    title = "Violin Plot: Weight by Gender",
    x = "Gender", y = "Weight (kg)"
  ) +
  scale_fill_manual(values = c("Female" = "#F8766D", "Male" = "#619CFF")) +
  theme_minimal()

# Pakker
library(ggplot2)
library(dplyr)

# Plot 1: Snout Length by Gender
ggplot(T6.19, aes(x = Gender, y = `Snout Vent Length (cm)`, fill = Gender)) +
  geom_violin(trim = FALSE, alpha = 0.6) +
  geom_boxplot(width = 0.1, fill = "white", outlier.shape = NA) +
  labs(
    title = "Violin Plot: Snout Length by Gender",
    x = "Gender",
    y = "Snout Length (cm)"
  ) +
  scale_fill_manual(values = c("Female" = "#F8766D", "Male" = "#619CFF")) +
  theme_minimal()

# Plot 2: Weight by Gender
ggplot(T6.19, aes(x = Gender, y = `Weight (kg)`, fill = Gender)) +
  geom_violin(trim = FALSE, alpha = 0.6) +
  geom_boxplot(width = 0.1, fill = "white", outlier.shape = NA) +
  labs(
    title = "Violin Plot: Weight by Gender",
    x = "Gender",
    y = "Weight (kg)"
  ) +
  scale_fill_manual(values = c("Female" = "#F8766D", "Male" = "#619CFF")) +
  theme_minimal()

  • Bruges til at supplere konklusionen fra multivariat test ————————————————————–

2.B Geometrisk forståelse af kovarians og korrelation (frivillig sektion)

# Udvælg de første 3 observationer (uanset køn)
X <- as.matrix(T6.19[1:3, 1:2])  # 3x2 matrix (Snout, Weight)

# Beregn middelvektor (kolonnegennemsnit)
x_bar <- colMeans(X)

# Lav 1-vektor
Ones <- rep(1, nrow(X))

# Beregn afvigelsesmatrix: Y = X - 1 * x_bar^T
Y <- X - Ones %*% t(x_bar)

# Udtræk kolonnerne som vektorer
Y1 <- Y[,1]  # Snout
Y2 <- Y[,2]  # Weight

# Beregn længder
Ly1 <- sqrt(sum(Y1^2))
Ly2 <- sqrt(sum(Y2^2))

# Beregn cosinus og vinkel
costh <- sum(Y1 * Y2) / (Ly1 * Ly2)
ang <- acos(costh) * 180 / pi

# Varians, kovarians og korrelation
s11 <- sum(Y1^2) / nrow(X)
s22 <- sum(Y2^2) / nrow(X)
s12 <- sum(Y1 * Y2) / nrow(X)
r12 <- s12 / sqrt(s11 * s22)

# Print resultater
cat("Afvigelsesvektorer for Anaconda-data (første 3 rækker):\n\n")
## Afvigelsesvektorer for Anaconda-data (første 3 rækker):
cat("Længde af afvigelsesvektor Snout:", Ly1, "\n")
## Længde af afvigelsesvektor Snout: 155.7996
cat("Længde af afvigelsesvektor Weight:", Ly2, "\n")
## Længde af afvigelsesvektor Weight: 50.37466
cat("Cosinus til vinkel mellem vektorer:", costh, "\n")
## Cosinus til vinkel mellem vektorer: 0.9957646
cat("Vinkel mellem vektorer (grader):", ang, "\n\n")
## Vinkel mellem vektorer (grader): 5.275185
cat("Varians af Snout:", s11, "\n")
## Varians af Snout: 8091.176
cat("Varians af Weight:", s22, "\n")
## Varians af Weight: 845.8689
cat("Kovarians:", s12, "\n")
## Kovarians: 2605.038
cat("Korrelationskoefficient:", r12, "\n")
## Korrelationskoefficient: 0.9957646
# Udtræk de første 3 observationer (Snout og Weight)
X <- as.matrix(T6.19[1:3, 1:2])
x_bar <- colMeans(X)
Ones <- rep(1, 3)

# Afvigelsesmatrix: Y = X - 1 * x_bar^T
Y <- X - Ones %*% t(x_bar)

# Udtræk kolonner
Y1 <- Y[,1]  # Snout
Y2 <- Y[,2]  # Weight

# Plot setup
plot(NA, xlim = range(c(0, Y1)), ylim = range(c(0, Y2)),
     xlab = "Snout deviation", ylab = "Weight deviation",
     main = "Afvigelsesvektorer fra origo (2D)", asp = 1)
grid()

# Tegn vektorer fra origo (0,0)
arrows(0, 0, Y1[1], Y2[1], col = "red", lwd = 2)
arrows(0, 0, Y1[2], Y2[2], col = "blue", lwd = 2)
arrows(0, 0, Y1[3], Y2[3], col = "darkgreen", lwd = 2)

# Tilføj labels
text(Y1[1], Y2[1], labels = "Y1", pos = 4, col = "red")
text(Y1[2], Y2[2], labels = "Y2", pos = 4, col = "blue")
text(Y1[3], Y2[3], labels = "Y3", pos = 4, col = "darkgreen")

3. Dekomponering af snout og vægt

# Udtræk de første 5 observationer for hver køn
T6.19_5females <- T6.19[T6.19$Gender == "Female", ][1:5, 1:2]
T6.19_5males   <- T6.19[T6.19$Gender == "Male",   ][1:5, 1:2]

# Kombiner i én matrix (2 rækker, 5 kolonner pr variabel)
snout_matrix <- rbind(
  T6.19_5males$`Snout Vent Length (cm)`,
  T6.19_5females$`Snout Vent Length (cm)`
)
rownames(snout_matrix) <- c("Male", "Female")

weight_matrix <- rbind(
  T6.19_5males$`Weight (kg)`,
  T6.19_5females$`Weight (kg)`
)
rownames(weight_matrix) <- c("Male", "Female")

# Midler
snout_mean <- matrix(rep(mean(snout_matrix), 10), nrow = 2, byrow = TRUE)
rownames(snout_mean) <- c("Male", "Female")

weight_mean <- matrix(rep(mean(weight_matrix), 10), nrow = 2, byrow = TRUE)
rownames(weight_mean) <- c("Male", "Female")

# Treatment effect
snout_group_means <- matrix(c(rep(mean(snout_matrix[1, ]), 5),
                              rep(mean(snout_matrix[2, ]), 5)),
                            nrow = 2, byrow = TRUE)
rownames(snout_group_means) <- c("Male", "Female")
snout_treatment <- snout_group_means - snout_mean
rownames(snout_treatment) <- c("Male", "Female")

weight_group_means <- matrix(c(rep(mean(weight_matrix[1, ]), 5),
                               rep(mean(weight_matrix[2, ]), 5)),
                             nrow = 2, byrow = TRUE)
rownames(weight_group_means) <- c("Male", "Female")
weight_treatment <- weight_group_means - weight_mean
rownames(weight_treatment) <- c("Male", "Female")

# Residual
snout_residual <- snout_matrix - snout_group_means
rownames(snout_residual) <- c("Male", "Female")

weight_residual <- weight_matrix - weight_group_means
rownames(weight_residual) <- c("Male", "Female")

# Vis resultat for Snout
cat("Snout Vent Length:\n\n")
## Snout Vent Length:
cat("Original data:\n"); print(snout_matrix)
## Original data:
##         [,1]  [,2]  [,3]  [,4] [,5]
## Male   176.7 259.5 258.0 229.8  233
## Female 271.0 477.0 306.3 365.3  466
cat("\nMean component:\n"); print(snout_mean)
## 
## Mean component:
##          [,1]   [,2]   [,3]   [,4]   [,5]
## Male   304.26 304.26 304.26 304.26 304.26
## Female 304.26 304.26 304.26 304.26 304.26
cat("\nTreatment effect:\n"); print(snout_treatment)
## 
## Treatment effect:
##          [,1]   [,2]   [,3]   [,4]   [,5]
## Male   -72.86 -72.86 -72.86 -72.86 -72.86
## Female  72.86  72.86  72.86  72.86  72.86
cat("\nResidual:\n"); print(snout_residual)
## 
## Residual:
##           [,1]  [,2]   [,3]   [,4]  [,5]
## Male    -54.70 28.10  26.60  -1.60  1.60
## Female -106.12 99.88 -70.82 -11.82 88.88
# Vis resultat for Weight
cat("\n\nWeight:\n\n")
## 
## 
## Weight:
cat("Original data:\n"); print(weight_matrix)
## Original data:
##        [,1]  [,2]  [,3] [,4]  [,5]
## Male    3.0  9.75 10.07  7.5  6.25
## Female 18.5 82.50 23.40 33.5 69.00
cat("\nMean component:\n"); print(weight_mean)
## 
## Mean component:
##          [,1]   [,2]   [,3]   [,4]   [,5]
## Male   26.347 26.347 26.347 26.347 26.347
## Female 26.347 26.347 26.347 26.347 26.347
cat("\nTreatment effect:\n"); print(weight_treatment)
## 
## Treatment effect:
##           [,1]    [,2]    [,3]    [,4]    [,5]
## Male   -19.033 -19.033 -19.033 -19.033 -19.033
## Female  19.033  19.033  19.033  19.033  19.033
cat("\nResidual:\n"); print(weight_residual)
## 
## Residual:
##           [,1]   [,2]    [,3]    [,4]   [,5]
## Male    -4.314  2.436   2.756   0.186 -1.064
## Female -26.880 37.120 -21.980 -11.880 23.620
# Som plot

library(tidyr)
library(dplyr)
library(ggplot2)

# ID og køn
ID <- paste0("Obs", 1:10)
Gender <- c(rep("Male", 5), rep("Female", 5))

# Rækkefølgebevarende faktor
ID <- factor(ID, levels = paste0("Obs", 1:10))

# ----- Snout -----
df_snout <- data.frame(
  ID = ID,
  Gender = Gender,
  Mean = as.vector(snout_mean),
  Treatment = as.vector(snout_treatment),
  Residual = as.vector(snout_residual)
) %>%
  pivot_longer(cols = c("Mean", "Treatment", "Residual"),
               names_to = "Component", values_to = "Value")

# ----- Weight -----
df_weight <- data.frame(
  ID = ID,
  Gender = Gender,
  Mean = as.vector(weight_mean),
  Treatment = as.vector(weight_treatment),
  Residual = as.vector(weight_residual)
) %>%
  pivot_longer(cols = c("Mean", "Treatment", "Residual"),
               names_to = "Component", values_to = "Value")

# ----- Plot Snout -----
ggplot(df_snout, aes(x = ID, y = Value, fill = Component)) +
  geom_bar(stat = "identity") +
  ggtitle("Dekomponering af Snout (Length)") +
  ylab("cm") +
  theme_minimal()

# ----- Plot Weight -----
ggplot(df_weight, aes(x = ID, y = Value, fill = Component)) +
  geom_bar(stat = "identity") +
  ggtitle("Dekomponering af Weight") +
  ylab("kg") +
  theme_minimal()

library(tidyr)
library(dplyr)
library(ggplot2)

# ID og køn
ID <- paste0("Obs", 1:10)
Gender <- c(rep("Male", 5), rep("Female", 5))
ID <- factor(ID, levels = paste0("Obs", 1:10))  # Rækkefølge

# ----- Snout -----
df_snout <- data.frame(
  ID = ID,
  Gender = Gender,
  Mean = as.vector(snout_mean),
  Treatment = as.vector(snout_treatment),
  Residual = as.vector(snout_residual)
) %>%
  pivot_longer(cols = c("Mean", "Treatment", "Residual"),
               names_to = "Component", values_to = "Value") %>%
  mutate(Component = factor(Component, levels = c("Mean", "Treatment", "Residual")))

# ----- Weight -----
df_weight <- data.frame(
  ID = ID,
  Gender = Gender,
  Mean = as.vector(weight_mean),
  Treatment = as.vector(weight_treatment),
  Residual = as.vector(weight_residual)
) %>%
  pivot_longer(cols = c("Mean", "Treatment", "Residual"),
               names_to = "Component", values_to = "Value") %>%
  mutate(Component = factor(Component, levels = c("Mean", "Treatment", "Residual")))

# ----- Plot Snout -----
ggplot(df_snout, aes(x = ID, y = Value, fill = Component)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = round(Value, 1)),
            position = position_stack(vjust = 0.5),
            size = 3, color = "black") +
  scale_fill_manual(values = c("Mean" = "#F8766D", "Treatment" = "#619CFF", "Residual" = "#00BA38")) +
  ggtitle("Dekomponering af Snout (Length)") +
  ylab("cm") +
  theme_minimal()

# ----- Plot Weight -----
ggplot(df_weight, aes(x = ID, y = Value, fill = Component)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = round(Value, 1)),
            position = position_stack(vjust = 0.5),
            size = 3, color = "black") +
  scale_fill_manual(values = c("Mean" = "#F8766D", "Treatment" = "#619CFF", "Residual" = "#00BA38")) +
  ggtitle("Dekomponering af Weight") +
  ylab("kg") +
  theme_minimal()

# Tilføj Gender label
gender <- c(rep("Female", 5), rep("Male", 5))

# Dataframes til snout
df_snout <- data.frame(
  ID = factor(paste0("Obs", 1:10), levels = paste0("Obs", 1:10)),
  Gender = gender,
  Mean = as.vector(snout_mean),
  Treatment = as.vector(snout_treatment),
  Residual = as.vector(snout_residual)
) %>%
  pivot_longer(cols = c("Mean", "Treatment", "Residual"),
               names_to = "Component", values_to = "Value")

# Dataframes til weight
df_weight <- data.frame(
  ID = factor(paste0("Obs", 1:10), levels = paste0("Obs", 1:10)),
  Gender = gender,
  Mean = as.vector(weight_mean),
  Treatment = as.vector(weight_treatment),
  Residual = as.vector(weight_residual)
) %>%
  pivot_longer(cols = c("Mean", "Treatment", "Residual"),
               names_to = "Component", values_to = "Value")

# ----- Plot Snout (med køn i farve) -----
ggplot(df_snout, aes(x = ID, y = Value, fill = Component)) +
  geom_bar(stat = "identity") +
  facet_wrap(~Gender, scales = "free_x") +
  ggtitle("Dekomponering af Snout (Length) opdelt efter køn") +
  ylab("cm") +
  theme_minimal()

# ----- Plot Weight (med køn i farve) -----
ggplot(df_weight, aes(x = ID, y = Value, fill = Component)) +
  geom_bar(stat = "identity") +
  facet_wrap(~Gender, scales = "free_x") +
  ggtitle("Dekomponering af Weight opdelt efter køn") +
  ylab("kg") +
  theme_minimal()

library(tidyr)
library(dplyr)
library(ggplot2)

# Tilføj Gender label
gender <- c(rep("Female", 5), rep("Male", 5))
ID <- factor(paste0("Obs", 1:10), levels = paste0("Obs", 1:10))

# ---- SNOUT ----
df_snout <- data.frame(
  ID = ID,
  Gender = gender,
  Mean = as.vector(snout_mean),
  Treatment = as.vector(snout_treatment),
  Residual = as.vector(snout_residual)
) %>%
  pivot_longer(cols = c("Mean", "Treatment", "Residual"),
               names_to = "Component", values_to = "Value") %>%
  mutate(Component = factor(Component, levels = c("Mean", "Treatment", "Residual")))

# ---- WEIGHT ----
df_weight <- data.frame(
  ID = ID,
  Gender = gender,
  Mean = as.vector(weight_mean),
  Treatment = as.vector(weight_treatment),
  Residual = as.vector(weight_residual)
) %>%
  pivot_longer(cols = c("Mean", "Treatment", "Residual"),
               names_to = "Component", values_to = "Value") %>%
  mutate(Component = factor(Component, levels = c("Mean", "Treatment", "Residual")))

# ---- Farver ----
custom_colors <- c("Mean" = "#F8766D", "Treatment" = "#619CFF", "Residual" = "#00BA38")

# ---- Plot SNOUT ----
ggplot(df_snout, aes(x = ID, y = Value, fill = Component)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = round(Value, 1)),
            position = position_stack(vjust = 0.5),
            size = 3, color = "black") +
  facet_wrap(~Gender, scales = "free_x") +
  scale_fill_manual(values = custom_colors) +
  ggtitle("Dekomponering af Snout (Length) opdelt efter køn") +
  ylab("cm") +
  theme_minimal()

# ---- Plot WEIGHT ----
ggplot(df_weight, aes(x = ID, y = Value, fill = Component)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = round(Value, 1)),
            position = position_stack(vjust = 0.5),
            size = 3, color = "black") +
  facet_wrap(~Gender, scales = "free_x") +
  scale_fill_manual(values = custom_colors) +
  ggtitle("Dekomponering af Weight opdelt efter køn") +
  ylab("kg") +
  theme_minimal()

4. Multivariat inferens

(Hotelling’s T² test: H₀: μ_female = μ_male)

4.1 Hotelling’s T² test * Hotelling’s T² test for forskel i middelvektor:

  • \(H_0: \mu_{\text{female}} = \mu_{\text{male}}\)
  • \(H_A: \mu_{\text{female}} \ne \mu_{\text{male}}\)
  • Hvis signifikant → fortsæt med simultane CI’er
# Udtræk data for hver gruppe (kun numeriske kolonner)
Xf <- T6.19[T6.19$Gender == "Female", 1:2]  # Snout og Weight
Xm <- T6.19[T6.19$Gender == "Male", 1:2]

# Gruppestørrelser
n1 <- nrow(Xf)
n2 <- nrow(Xm)
p <- ncol(Xf)  # Antal variable

# Gruppemiddelværdier
xbar1 <- colMeans(Xf)
xbar2 <- colMeans(Xm)

# Pooled kovariansmatrix
S1 <- cov(Xf)
S2 <- cov(Xm)
Spooled <- ((n1 - 1)*S1 + (n2 - 1)*S2)/(n1 + n2 - 2)

# Hotelling’s T² teststørrelse
T2 <- (n1 * n2) / (n1 + n2) * t(xbar1 - xbar2) %*% solve(Spooled) %*% (xbar1 - xbar2)

# Kritisk værdi og p-værdi
F_crit <- ((n1 + n2 - 2) * p / (n1 + n2 - p - 1)) * qf(0.95, p, n1 + n2 - p - 1)
F_val <- ((n1 + n2 - p - 1) * T2) / ((n1 + n2 - 2) * p)
p_val <- 1 - pf(F_val, p, n1 + n2 - p - 1)

# Udskriv resultat
cat("Hotelling’s T² test:\n")
## Hotelling’s T² test:
cat("T² =", round(T2, 4), "\n")
## T² = 76.9153
cat("F =", round(F_val, 4), "\n")
## F = 37.7455
cat("p-værdi =", format.pval(p_val, digits = 5), "\n")
## p-værdi = 0.000000000064296

Hotelling’s T² test undersøger, om der er en samlet forskel i middelvektorerne for Snout Length og Weight mellem hunner og hanner.

Da p-værdien er ekstremt lav (p ≪ 0.05), forkaster vi nulhypotesen H₀: μ_female = μ_male.

Konklusion: Der er en signifikant multivariat forskel på længde og vægt mellem de to køn.

4.2 Generaliseret stikprøvevarians (|S|)

Generaliseret stikprøvevarians (|S|) for poolet kovariansmatrix

Total Spredning i Flere Dimensioner (|S|)

Multivariat Spredning: Determinant af Spooled

kovariansstruktur: Generaliseret Varians

Determinant af Spooled Matrix – Mål for Multivariat Spredning

|S| som Mål for Volumen af Dataspredning

(relevant for Hotelling’s T² og Mahalanobis-afstand)

# 1. Generaliseret stikprøvevarians for pooled matrix
Xf <- T6.19[T6.19$Gender == "Female", 1:2]
Xm <- T6.19[T6.19$Gender == "Male", 1:2]
n1 <- nrow(Xf); n2 <- nrow(Xm)
S1 <- cov(Xf); S2 <- cov(Xm)
Spooled <- ((n1 - 1)*S1 + (n2 - 1)*S2)/(n1 + n2 - 2)

# Determinant = generaliseret varians
gen_var <- det(Spooled)
cat("Generaliseret stikprøvevarians (|S|):", round(gen_var, 4))
## Generaliseret stikprøvevarians (|S|): 86170.38

Fortolkning

Ekstra fortolkning til posteren:

cat("\n|S| udtrykker den samlede multivariate variation og volumen i datasættet.\n",
    "En høj værdi indikerer stor spredning i både længde og vægt.\n")
## 
## |S| udtrykker den samlede multivariate variation og volumen i datasættet.
##  En høj værdi indikerer stor spredning i både længde og vægt.

|S| udtrykker den samlede multivariate variation og volumen i datasættet. En høj værdi indikerer stor spredning i både længde og vægt.

cat("Generaliseringen |S| = |D| ⋅ |R| viser, at varians og korrelation sammen bestemmer spredningens volumen.\n")
## Generaliseringen |S| = |D| ⋅ |R| viser, at varians og korrelation sammen bestemmer spredningens volumen.
cat("Bemærk: Determinanten af Spooled-matrixen fungerer som en generaliseret stikprøvevarians.\n",
    "En lav værdi tyder på høj korrelation mellem variablerne, mens en høj værdi tyder på større uafhængighed.\n")
## Bemærk: Determinanten af Spooled-matrixen fungerer som en generaliseret stikprøvevarians.
##  En lav værdi tyder på høj korrelation mellem variablerne, mens en høj værdi tyder på større uafhængighed.

|S| måler samlet spredning. Lav værdi → høj korrelation. Høj værdi → større uafhængighed.

4.3 Visualisering: Hotelling’s T²-ellipse

#Grafisk med ellipse pakken

library(ellipse)
## 
## Attaching package: 'ellipse'
## The following object is masked from 'package:car':
## 
##     ellipse
## The following object is masked from 'package:graphics':
## 
##     pairs
# Brug Spooled og center i forskellen mellem køn
Xf <- T6.19[T6.19$Gender == "Female", 1:2]
Xm <- T6.19[T6.19$Gender == "Male", 1:2]
xbar1 <- colMeans(Xf)
xbar2 <- colMeans(Xm)
diff_mean <- xbar1 - xbar2

n1 <- nrow(Xf); n2 <- nrow(Xm)
Spooled <- ((n1 - 1)*cov(Xf) + (n2 - 1)*cov(Xm))/(n1 + n2 - 2)
n <- (n1 * n2) / (n1 + n2)
alpha <- 0.05
p <- 2
T2_crit <- ((n1 + n2 - 2) * p / (n1 + n2 - p - 1)) * qf(1 - alpha, p, n1 + n2 - p - 1)

# Plot
plot(ellipse(Spooled / n, centre = diff_mean, level = 1 - alpha,
             t = sqrt(T2_crit)), type = "l", col = "darkblue",
     xlab = "Snout Vent Length (cm)", ylab = "Weight (kg)",
     main = "Hotelling’s T² Ellipse (Volumen ≈ |S|)")
points(0, 0, pch = 4, col = "red")
text(0, 0, "0", pos = 1, cex = 0.8)

# Grafisk uden ellipse pakken
# Parametre
xbar1 <- colMeans(Xf)
xbar2 <- colMeans(Xm)
diff_mean <- xbar1 - xbar2
n <- (n1 * n2) / (n1 + n2)
alpha <- 0.05
p <- 2
T2_crit <- ((n1 + n2 - 2) * p / (n1 + n2 - p - 1)) * qf(1 - alpha, p, n1 + n2 - p - 1)

# Generér ellipse uden ellipse-pakke
theta <- seq(0, 2*pi, length.out = 200)
circle <- cbind(cos(theta), sin(theta))
A <- chol(Spooled / n * T2_crit)
ellipse_coords <- t(diff_mean + t(circle %*% A))

# Plot
plot(ellipse_coords, type = "l", lwd = 2, col = "black",
     xlab = "Snout Vent Length (cm)", ylab = "Weight (kg)",
     main = "Hotelling’s T² Ellipse (uden ellipse-pakke)")
points(0, 0, col = "blue", pch = 19)
text(0, 0, "0", pos = 1)
points(diff_mean[1], diff_mean[2], pch = 19, col = "darkblue")
#text(diff_mean[1], diff_mean[2], "M1 - M2", pos = 4, col = "darkblue")
text(diff_mean[1], diff_mean[2], expression(mu[1] - mu[2]), pos = 4, col = "darkblue")

🔍 Punktet “M1 - M2” repræsenterer forskellen mellem de to gruppers stikprøvegennemsnit. Hvis dette punkt falder uden for T²-ellipsen, tyder det på en signifikant forskel i middelværdierne.

✅ Kort opsummering til posteren:

Den sorte ellipse repræsenterer området for ikke-signifikante differenser i middelvektorer mellem de to køn.

Da det blå punkt (μ₁ − μ₂) falder uden for ellipsen, forkastes H₀ → der er en signifikant forskel.

4.4 Mannova (Simpel med wilks)

# Multivariat variansanalyse (MANOVA)
model_manova <- manova(cbind(`Snout Vent Length (cm)`, `Weight (kg)`) ~ Gender, data = T6.19)

# Samlet test med summary
summary(model_manova, test = "Wilks")
##           Df   Wilks approx F num Df den Df          Pr(>F)    
## Gender     1 0.41248   37.745      2     53 0.0000000000643 ***
## Residuals 54                                                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Kort forklaring til posteren (valgfrit)

Der blev udført en multivariat variansanalyse (MANOVA) for at teste, om der er en samlet forskel i længde og vægt mellem hanner og hunner.
Wilks’ lambda-test viste en meget signifikant forskel (p < 2.2e-16), hvilket bekræfter resultaterne fra de univariate t-tests.

4.4b Mannova med Bartletts korrektion

library(dplyr)

# Data: responsvariabler og grupper
Y <- T6.19[, 1:2]  # Snout & Weight
group <- T6.19$Gender
n1 <- sum(group == "Female")
n2 <- sum(group == "Male")
g <- 2
p <- ncol(Y)
n_total <- nrow(Y)

# Fit MANOVA-model
fit <- manova(as.matrix(Y) ~ group)

# Kort metode – direkte Wilks-test
summary(fit, test = "Wilks")
##           Df   Wilks approx F num Df den Df          Pr(>F)    
## group      1 0.41248   37.745      2     53 0.0000000000643 ***
## Residuals 54                                                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Lang metode – Udtræk SSCP-matricer
B <- summary(fit)$SS$group
W <- summary(fit)$SS$Residuals
T_mat <- B + W

# Frihedsgrader
df1 <- p * (g - 1)          # for behandling
df2 <- p * (n_total - g)    # for residualer

# Wilks' lambda
lambda <- det(W) / det(T_mat)

# Chi-i-anden Bartletts korrektion (Resultat 6.8 / 6.43)
m <- (n1 + n2 - 2) - (p + g) / 2
chi2_stat <- -m * log(lambda)
chi2_crit <- qchisq(0.95, df1)
pval_chi2 <- pchisq(chi2_stat, df1, lower.tail = FALSE)

# Udskriv resultater
cat("Wilks’ lambda: ", round(lambda, 4), "\n")
## Wilks’ lambda:  0.4125
cat("Chi²-test (Bartletts korrektion):\n")
## Chi²-test (Bartletts korrektion):
cat("Teststatistik:", round(chi2_stat, 3), "\n")
## Teststatistik: 46.049
cat("Kritisk værdi (χ², 95%, df =", df1, "):", round(chi2_crit, 3), "\n")
## Kritisk værdi (χ², 95%, df = 2 ): 5.991
cat("p-værdi:", signif(pval_chi2, 4), "\n")
## p-værdi: 0.0000000001001
cat("Konklusion: ", ifelse(pval_chi2 < 0.05, "Forkast H₀ – Køn påvirker mindst én af variablerne.", "Forkast ikke H₀ – Ingen forskel."), "\n")
## Konklusion:  Forkast H₀ – Køn påvirker mindst én af variablerne.

📌 MANOVA-resultat (ud fra Bartletts korrektion)

Wilks’ λ = 0.4125

χ² = 46.05, kritisk χ² = 5.99, df = 2

p < 0.0000000001

✅ Konklusion: Forkast H₀ — køn har signifikant indflydelse på mindst én af variablerne (Snout Length eller Weight).

SAmlet konklusion på MANNVA:

📌 Kort fortolkning til posteren (tekst)

MANOVA: Multivariat variansanalyse

Undersøger om middelvektorerne for Snout Length og Weight er ens mellem hunner og hanner:

    H₀: μ<sub>female</sub> = μ<sub>male</sub>

    H₁: μ<sub>female</sub> ≠ μ<sub>male</sub>

Resultat:

    Wilks’ λ = 0.4125

    F(2, 53) = 37.745, p < 0.0000000001

    χ² = 46.05, kritisk χ² = 5.99, df = 2

    p < 0.0000000001

✅ Begge tests fører til samme konklusion:
Forkast H₀ — køn påvirker signifikant mindst én af variablerne (snout eller vægt eller begge).

5. Test for kovariansmatricer (Box’s M test)

# Box’s M test for lighed i kovariansmatricer

# Forudsæt: T6.19 har kolonner: Snout Vent Length (cm), Weight (kg), Gender
library(dplyr)

# Del op efter køn
Xf <- T6.19 %>%
  filter(Gender == "Female") %>%
  dplyr::select(`Snout Vent Length (cm)`, `Weight (kg)`)

Xm <- T6.19 %>%
  filter(Gender == "Male") %>%
  dplyr::select(`Snout Vent Length (cm)`, `Weight (kg)`)


n1 <- nrow(Xf)
n2 <- nrow(Xm)
p <- ncol(Xf)
g <- 2
n_total <- n1 + n2

# Kovariansmatricer
S1 <- cov(Xf)
S2 <- cov(Xm)

# Determinanter
detS1 <- det(S1)
detS2 <- det(S2)

# Spooled
Spooled <- ((n1 - 1) * S1 + (n2 - 1) * S2) / (n_total - g)
detSpooled <- det(Spooled)

# Box's M test statistik
M <- (n1 - 1) * log(detSpooled / detS1) + (n2 - 1) * log(detSpooled / detS2)

# Korrektionsfaktor
u <- ((1 / (n1 - 1)) + (1 / (n2 - 1)) - (1 / (n_total - g))) *
     (2 * p^2 + 3 * p - 1) / (6 * (p + 1) * (g - 1))

C <- (1 - u) * M
df <- p * (p + 1) * (g - 1) / 2
pval <- pchisq(C, df, lower.tail = FALSE)

# Resultat
cat("Box’s M test:\n")
## Box’s M test:
cat("Teststatistik:", round(C, 3), "\n")
## Teststatistik: 100.325
cat("Frihedsgrader:", df, "\n")
## Frihedsgrader: 3
cat("p-værdi:", signif(pval, 4), "\n")
## p-værdi: 0.000000000000000000001323
cat("Konklusion: ", ifelse(pval < 0.05,
                          "Forkast H₀ – kovariansmatricerne er forskellige.",
                          "Forkast ikke H₀ – antag ens kovariansmatricer."), "\n")
## Konklusion:  Forkast H₀ – kovariansmatricerne er forskellige.

📌 2. Fortolkning til posteren (Box’s M-test)

Box’s M-test undersøger, om varians-/kovariansstrukturen er ens mellem grupperne — her: hunner og hanner. Hypoteser:

H₀: Kovariansmatricerne for de to køn er ens

H₁: Mindst én gruppe har en anden kovariansstruktur

Resultat:

Teststatistik: 100.325

Frihedsgrader: 3

p-værdi: < 0.0000000000000000000014

✅ Konklusion: Da p-værdien er langt under 0.05, forkastes H₀. Der er signifikant forskel i kovariansmatricerne — dvs. variabiliteten i længde og vægt varierer mellem hunner og hanner.

6. Simultane konfidensintervaller

(Ellipse + Bonferroni CI)

library(MASS)

# Del data
Xf <- T6.19[T6.19$Gender == "Female", 1:2]
Xm <- T6.19[T6.19$Gender == "Male", 1:2]

# Parametre
xbar1 <- colMeans(Xf)
xbar2 <- colMeans(Xm)
diff_mean <- xbar1 - xbar2
n1 <- nrow(Xf); n2 <- nrow(Xm)
p <- ncol(Xf)
n <- (n1 * n2) / (n1 + n2)
Spooled <- ((n1 - 1)*cov(Xf) + (n2 - 1)*cov(Xm)) / (n1 + n2 - 2)

# Kritiske værdier
alpha <- 0.05
t_crit <- qt(1 - alpha/(2*p), df = n1 + n2 - 2)  # Bonferroni
T2_crit <- ((n1 + n2 - 2)*p / (n1 + n2 - p - 1)) * qf(1 - alpha, p, n1 + n2 - p - 1)  # T²

# Bonferroni-interval (rektangel)
margin <- t_crit * sqrt(diag(Spooled) / n)
rect_coords <- rbind(
  diff_mean + c(-margin[1], -margin[2]),
  diff_mean + c( margin[1], -margin[2]),
  diff_mean + c( margin[1],  margin[2]),
  diff_mean + c(-margin[1],  margin[2]),
  diff_mean + c(-margin[1], -margin[2])
)

# T²-ellipse
theta <- seq(0, 2*pi, length.out = 200)
circle <- cbind(cos(theta), sin(theta))
A <- chol(Spooled / n * T2_crit)
ellipse_coords <- t(diff_mean + t(circle %*% A))

# Plot
plot(ellipse_coords, type = "l", col = "blue", lwd = 2,
     xlab = "Snout Vent Length (cm)", ylab = "Weight (kg)",
     main = "Konfidensellipse (T²) og Bonferroni-rektangel")
lines(rect_coords, col = "red", lty = 2, lwd = 2)
points(0, 0, pch = 4, col = "black")
text(0, 0, "0", pos = 1)
points(diff_mean[1], diff_mean[2], col = "blue", pch = 19)
text(diff_mean[1], diff_mean[2], expression(mu[1] - mu[2]), pos = 4, cex = 0.8)
legend("topright", legend = c("T²-ellipse", "Bonferroni-interval"),
       col = c("blue", "red"), lty = c(1, 2), lwd = 2)

🧠 2. Fortolkning – til posteren

📌 Konfidensellipse og -rektangel

Vi sammenligner forskellen mellem hunner og hanner ift. både længde og vægt.

To konfidensområder vises:

    🔵 Ellipse (Hotelling’s T²): Multivariat konfidensområde

    🔺 Rektangel (Bonferroni CI): To separate simultane konfidensintervaller

0,0-punktet repræsenterer ingen forskel mellem kønnene.

Da både ellipse og rektangel ikke indeholder 0,0 → vi har stærk evidens for forskel i mindst én variabel.

7. Konklusion og fortolkning

📌 Diskussion og konklusioner

Gennem både univariate og multivariate analyser fandt vi klare tegn på, at køn har signifikant indflydelse på både Snout Length og Weight hos anacondaerne.

T-tests viste signifikante forskelle i både længde og vægt mellem hunner og hanner (p ≪ 0.05).

Violinplots og boxplots understøttede disse resultater visuelt med tydelig adskillelse mellem kønnene.

Hotelling’s T²-test og MANOVA bekræftede, at forskellen i middelvektorer er signifikant (Wilks' λ = 0.4125, p < 0.0000000001).

Simultane konfidensintervaller (Bonferroni og T²-ellipse) viste begge, at nulpunktet (0,0) ikke ligger inde i intervallerne → dvs. vi kan med høj sikkerhed afvise, at der ikke er forskel.

Box’s M-test viste, at variansen og kovariansstrukturen ikke er ens mellem kønnene (p < 0.000000000001), hvilket understreger forskellene yderligere.

🔍 Overordnet set konkluderer vi:

Der er statistisk belæg for, at køn har en signifikant effekt på både længde og vægt hos anacondaerne, og denne effekt ses både i gennemsnit, variation og kovariansstruktur.

Disse resultater bør tages i betragtning i biologisk forskning, fx ved artsopdeling, vækststudier eller populationsanalyse.